function osl_test_script_group(testdatadir,testoutputdir,S)

global OSLDIR;

testoutputdir=[testoutputdir '_group'];
runcmd(['rm -rf ' testoutputdir]);
mkdir(testoutputdir);

testplotsdir=[testoutputdir '/plots'];
mkdir(testplotsdir);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% osl_example_group_sensorspace_oat

printprefix='group_sensorspace_oat';
printindex=1;

%%%
% load pre-run OAT analysis and run group-level stage for SENSOR space analysis

oatdir=[testdatadir '/faces_group_data_new/faces_group_sensor_norm.oat'];
oat = osl_load_oat(oatdir); 

oat.do_plots=0;
oat.to_do=[0 0 0 1];

oat = osl_run_oat(oat);

S2=[];
S2.oat=oat;
S2.stats_fname=oat.group_level.results_fnames;
S2.modality='MEGPLANAR';
S2.first_level_contrast=1:3;
S2.group_level_contrast=1;
S2.cfg.colorbar='yes';
S2.view_cope=0;
S2.cfg.interactive = 'no';

% calculate t-stat using contrast of absolute value of parameter estimates
[cfg, data]=osl_stats_multiplotER(S2);

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% osl_example_group_oat

printprefix='group_oat';
printindex=1;

oatdir=[testdatadir '/faces_group_data_new/faces_group_sensor_norm.oat'];

%testdir=[tilde '/Desktop/faces_group_data_new']; 
%oatdir=[testdir '/faces_group_norm.oat'];

% Load OAT analysis for which the first 2 stages (source recon and
% first-level GLM) have already been run.
% Note that the 1st level contrasts run are:
% S2.contrast{1}=[3 0 0 0]'; % motorbikes
% S2.contrast{2}=[0 1 1 1]'; % faces
% S2.contrast{3}=[-3 1 1 1]'; % faces-motorbikes

oat = osl_load_oat(oatdir,'first_level','sub_level','group_level'); 

%% run group-level stage

oat.group_level.time_range=[0 0.3];
oat.group_level.space_average=0;
oat.group_level.time_average=0;
oat.group_level.time_smooth_std=0; % secs
oat.group_level.spatial_smooth_fwhm=0; % mm
oat.group_level.mask_fname='';oat.group_level=rmfield(oat.group_level,'mask_fname');
oat.group_level.use_tstat=1;
oat.group_level.group_varcope_time_smooth_std=0;
oat.group_level.group_varcope_spatial_smooth_fwhm=100;
oat.group_level.name='group_level';

% run OAT
oat.to_do=[0 0 0 1]; % run group-level stage only

oat = osl_check_oat(oat);

oat = osl_run_oat(oat);

% OUTPUT GROUP'S NIFTII FILES

S2=[];
S2.oat=oat;
S2.stats_fname=oat.group_level.results_fnames;
S2.first_level_contrasts=[3]; % list of first level contrasts to output

[statsdir,times]=osl_save_nii_stats(S2);

% con=3;runcmd(['fslview ' OSLDIR '/std_masks/MNI152_T1_2mm_brain.nii.gz ' OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex ' statsdir '/tstat' num2str(con) '_gc1_2mm.nii.gz ' statsdir '/cope' num2str(con) '_gc1_2mm.nii.gz &']);

figure;
con=3;
map=ra([statsdir '/tstat' num2str(con) '_gc1_2mm']);
resamp_gridstep=2;
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
percfrom=99;percto=99.99;
low=percentile((x1),percfrom);high=percentile((x1),percto);
xslice=49;vol=18;
overlay_act(flipud(squeeze(map(xslice,:,:,vol))'), flipud(squeeze(bgmap(xslice,:,:))'),'red2yellow',0,[low high],[3000 8000]);
text(2,2,[num2str(low) ',' num2str(high)],'Color','w');
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;
percfrom=99.5;percto=99.99;
low=percentile((x1),percfrom);high=percentile((x1),percto);
zslice=33;vol=18;
overlay_act(flipud(squeeze(map(:,:,zslice,vol))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
text(2,2,[num2str(low) ',' num2str(high)],'Color','w');
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,20:2:70,vol));colorbar;title('Group ERF: tstat3');
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;
con=3;
resamp_gridstep=2;
map=ra([statsdir '/cope' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
percfrom=99.6;percto=100;
low=percentile((x1),percfrom);high=percentile((x1),percto);
zslice=33;vol=18;
overlay_act(flipud(squeeze(map(:,:,zslice,vol))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
text(2,2,[num2str(low) ',' num2str(high)],'Color','w');
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,20:2:70,vol));colorbar;title('Group ERF: cope3');
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;


%%%%%%%%%%%%%%%%%%
%% Using ROI at the output from the group level

% load OAT analysis for which the first 3 stages have already been run
oat = osl_load_oat(oatdir,'first_level','sub_level','group_level'); 

mni_coords=[20,-68,-6]; %RFFA

S2=[];
S2.oat=oat;
S2.stats_fname=oat.group_level.results_fnames;
%S2.mask_fname=osl_mnicoords2mnimask(mni_coords,1,'mask');
S2.mni_coords=mni_coords;
[stats,times,mni_coords_used]=osl_output_roi_stats(S2);

tim=times;

baseline_correct=0;
gcon=1;
figure;ho;
cons=[1:3];
cols={'r','g','b'};
for coni=1:length(cons),
    con=cons(coni);
    cope=permute(stats.cope(1,:,con,:,gcon),[2 1 3 4 5]);
    stdcope=permute(stats.stdcope(1,:,con,:,gcon),[2 1 3 4 5]);
    mn=0;if(baseline_correct),mn=mean(cope(tim>-0.1 & tim<0));end;
    errorbar(tim,cope-mn,stdcope,cols{coni});
end;
plot4paper('time (secs)','COPE'); 
title('Group average copes (ERFs)');
legend('motorbikes','faces','faces-motorbikes'); 

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

% tstat
baseline_correct=0;
gcon=1;
vox_coord=1;
figure;ho;
cons=[1:3];
cols={'r','g','b'};
for coni=1:length(cons),
    con=cons(coni);
    cope=permute(stats.cope(1,:,con,:,gcon),[2 1 3 4 5]);
    stdcope=permute(stats.stdcope(1,:,con,:,gcon),[2 1 3 4 5]);
    mn=0;if(baseline_correct),mn=mean(cope(tim>-0.1 & tim<0));end;
    plot(tim,(cope-mn)./stdcope,cols{coni});
end;
plot4paper('time (secs)','t-stat'); 
title('Group average t-stats (ERFs)');
legend('motorbikes','faces','faces-motorbikes'); 

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

% figure;plot(tim,squeeze(mean(tmp(vox_coord,:,:),2))./sqrt(squeeze(var(tmp(vox_coord,:,:),[],2)./(size(tmp,2)-1)))); title('tstat');

%%%%%%%%%%%%%%%%%%
%% Using ROI at the input to the group level

% load OAT analysis for which the first 3 stages have already been run
oat = osl_load_oat(oatdir,'first_level','sub_level','group_level'); 

% mask
oat.group_level.time_range=[-0.15 0.3];
oat.group_level.space_average=1;
oat.group_level.time_average=0;
oat.group_level.time_smooth_std=0; % secs
oat.group_level.group_varcope_time_smooth_std=0.2;
oat.group_level.group_varcope_spatial_smooth_fwhm=100;
oat.group_level.spatial_smooth_fwhm=0; % mm
oat.group_level.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];

%oat.group_level.mask_fname=['/Users/woolrich/vols_data/murphy_project/data/masks/Right_Amygdala'];

oat.group_level.name='roi_average_group_level';
oat.group_level.use_tstat=1;

oat = osl_check_oat(oat);

oat.to_do=[0 0 0 1];

% run OAT
oat = osl_run_oat(oat);

% do plots

gstats=osl_load_oat_results(oat,oat.group_level.results_fnames);

tim=gstats.times
cons=[1:3];

figure;
baseline_correct=0;
gcon=1;
subplot(1,2,1);ho;
cols={'r','g','b'};
for coni=1:length(cons),
    con=cons(coni);
    cope=(permute(mean(gstats.cope(:,:,con,:,gcon),1),[2 4 1 3 5]));
    stdcope=(permute(mean(gstats.stdcope(:,:,con,:,gcon),1),[2 4 1 3 5]));
    mn=0;if(baseline_correct),mn=mean(cope(tim>-0.1 & tim<0));end;
    errorbar(tim,cope-mn,stdcope,cols{coni});
end;
plot4paper('time (secs)','COPE'); 
title(['Group average copes (ERFs), use_tstat=' num2str(oat.group_level.use_tstat)]);
legend('motorbikes','faces','faces-motorbikes'); 
%exportfig(gcf,[murphy_plots_dir '/' Sin.resdir_postfix '_con' num2str(Sin.contrast4localisation) '_' Sin.mask_fname '/p_fan_cope'],'Color','rgb','Format','eps');
 
% tstat
baseline_correct=0;
gcon=1;
subplot(1,2,2);ho;
cols={'r','g','b'};
for coni=1:length(cons),
    con=cons(coni);
    cope=permute(mean(gstats.cope(1,:,con,:,gcon),1),[2 4 1 3 5]);
    stdcope=permute(mean(gstats.stdcope(1,:,con,:,gcon),1),[2 4 1 3 5]);
    mn=0;if(baseline_correct),mn=mean(cope(tim>-0.1 & tim<0));end;
    plot(tim,(cope-mn)./stdcope,cols{coni});
end;
plot4paper('time (secs)','t-stat'); 
title(['Group average t-stats (ERFs), use_tstat=' num2str(oat.group_level.use_tstat)]);
legend('motorbikes','faces','faces-motorbikes'); 

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

% individual subject plots:

baseline_correct=1;
con=3;
tmp=((gstats.lower_level_copes{con}));
mn=0;if(baseline_correct),mn=mean(squash(tmp(1,:,tim<0)));end;
figure;plot(tim,permute(tmp(1,:,:),[2,3,1])-mn);ho;
plot(tim,squeeze(median(tmp(1,:,:)-mn,2)),'LineWidth',3);
plot4paper('time (secs)','1st level cope'); 
title(['First level cope ' num2str(con) ', use tstat=' num2str(oat.group_level.use_tstat)]);

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

% figure;plot(tim,squeeze(mean(tmp(1,:,:),2))./sqrt(squeeze(var(tmp(1,:,:),[],2)./(size(tmp,2)-1)))); title('tstat');

%%%%%%%%%%%%%%%%%%%%%%%%
%% Produce a single volume average over a time window

% load OAT analysis for which the first 2 stages have already been run
oat = osl_load_oat(oatdir,'first_level','sub_level','group_level'); 

gstats=osl_load_oat_results(oat,oat.group_level.results_fnames);

S=[];
S.data=gstats;
S.resample_method='collapse';
S.reduce_time=1;
S.time_range=[0.140 0.150];
[data] = osl_reduce_data_to_visualize(S)

S2=[]; 
S2.oat=oat;
S2.stats=data;
S2.first_level_contrasts=[3]; % list of first level contrasts to output
S2.stats_dir=[oat.source_recon.dirname '/' data.fname '_time_collapsed_stats'];
[statsdir,times]=osl_save_nii_stats(S2);

%runcmd(['fslview ' OSLDIR '/std_masks/MNI152_T1_2mm_brain.nii.gz ' OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex ' statsdir '/tstat3_gc1_2mm.nii.gz ' statsdir '/cope3_gc1_2mm.nii.gz &']);


figure;
con=3;
map=ra([statsdir '/tstat' num2str(con) '_gc1_2mm']);
resamp_gridstep=2;bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
percfrom=95;percto=99.99;
low=percentile((x1),percfrom);high=percentile((x1),percto);
zslice=30;vol=1;
overlay_act(flipud(squeeze(map(:,:,zslice,vol))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
text(2,2,[num2str(low) ',' num2str(high)],'Color','w');
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%% Permutation stats on a single volume

S=[];
S.oat=oat;
S.timepoint=[nearest(gstats.times,0.144)];
S.cluster_stats_thresh=6;
S.cluster_stats_nperms=200;
S.first_level_copes_to_do=[3];
S.group_level_copes_to_do=[1];
S.group_varcope_spatial_smooth_fwhm=100;
S.write_cluster_script=0;

[ gstats, statsdir ] = osl_cluster_permutation_testing( S );

% View permutation stats
% runcmd(['fslview ' OSLDIR '/std_masks/MNI152_T1_2mm_brain.nii.gz ' OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex ' gstats.dir '/tstat3_gc1_2mm.nii.gz ' gstats.dir '/clustere_tstat3_gc1_2mm.nii.gz ' gstats.dir '/clustere_corrp_tstat3_gc1_2mm.nii.gz &']);

figure;
con=3;
map=ra([statsdir '/clustere_corrp_tstat3_gc1_2mm']);
resamp_gridstep=2;bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
low=0.95; high=1;
zslice=30;vol=1;
overlay_act(flipud(squeeze(map(:,:,zslice,vol))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
text(2,2,[num2str(low) ',' num2str(high)],'Color','w');
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

if(0),

    %% Permutation stats in 4D

    S=[];
    S.oat=oat;
    %S.timepoint=[nearest(gstats.times,0.144)];
    
    S.cluster_stats_thresh=6;
    S.cluster_stats_nperms=200;
    S.first_level_copes_to_do=[3];
    S.group_level_copes_to_do=[1];
    S.group_varcope_spatial_smooth_fwhm=100;
    S.write_cluster_script=0;
    S.matlab_exe_name='/Applications/MATLAB_R2012a.app/bin/matlab';
    [ gstats ] = osl_cluster_permutation_testing( S );

    % VIEW NIFTII RESULTS IN FSLVIEW
    S2=[];
    S2.oat=oat;
    S2.stats_fname=oat.group_level.results_fnames;
    S2.first_level_contrasts=[3]; % list of first level contrasts to output
    [statsdir,times]=osl_save_nii_stats(S2);

    %con=3;
    %runcmd(['fslview ' OSLDIR '/std_masks/MNI152_T1_2mm_brain.nii.gz ' OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex ' statsdir '/tstat' num2str(con) '_gc1_2mm.nii.gz ' statsdir '/tstat' num2str(con) '_gc1_clust4d_corrp_2mm.nii.gz ' statsdir '/tstat' num2str(con) '_gc1_clust4d_2mm.nii.gz &']);
   

    figure;
    con=3;
    map=ra([statsdir '/tstat' num2str(con) '_gc1_2mm']);
    resamp_gridstep=2;bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
    x1=squash(map,abs(map));
    percfrom=99.5;percto=99.9;
    low=percentile((x1),percfrom);high=percentile((x1),percto);
    zslice=30;vol=61;
    overlay_act(flipud(squeeze(map(:,:,zslice,vol))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
    text(2,2,[num2str(low) ',' num2str(high)],'Color','w');
    print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

end;

%%%%%%%%%%%%%%%%%%%
%% ROI Time-Freq analysis

oat.source_recon.dirname=[testdir '/faces_group_norm_roi.oat'];
oat=osl_load_oat(oat.source_recon.dirname,'first_level_tf','sub_level','group_level');

oat.to_do=[0 0 0 1];

oat.group_level.time_range=oat.first_level.time_range;
oat.group_level.space_average=1;
oat.group_level.time_average=0;
oat.group_level.time_smooth_std=0; % secs
oat.group_level.spatial_smooth_fwhm=0; % mm
oat.group_level.use_tstat=0;
oat.group_level.name='group_level';

oat=osl_check_oat(oat);

oat=osl_run_oat(oat);

%% Plot the group results

% load GLM result
stats=osl_load_oat_results(oat,oat.group_level.results_fnames);

% Plot the parameter estimates for our closest beamformed voxel over time, for 
% our two contrasts of interest:
figure;
con=3; imagesc(stats.times, stats.frequencies, squeeze((stats.cope(1,:,con,:)))');axis xy;
ylabel('frequency (Hz)'); xlabel('time (s)'); colorbar; title(['cope' num2str(con)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;
imagesc(stats.times, stats.frequencies, squeeze(stats.cope(1,:,con,:)./stats.stdcope(1,:,con,:))',[-4 4]);axis xy;
ylabel('frequency (Hz)'); xlabel('time (s)'); colorbar; title(['tstat' num2str(con)]);

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;


% Plot the time course of a single freq bin:
figure;
freqbinfrom=nearest(stats.frequencies,13);
freqbinto=nearest(stats.frequencies,15);

freqbin=freqbinfrom:freqbinto

cons=1:3;cols={'r','g','b'};
for c=1:length(cons),
    con=cons(c); plot(stats.times, squeeze(mean(stats.cope(1,:,con,freqbin)./stats.stdcope(1,:,con,freqbin),4)),cols{c},'LineWidth',2); hold on;
end;
legend('Motorbikes','Faces','Faces-Motorbikes');
xlabel('time (s)'); ylabel('1-tailed t-stat'); title([num2str(stats.frequencies(freqbin)) ' Hz']);

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;


%% Do cluster permutation stats

% load GLM result
stats=osl_load_oat_results(oat,oat.group_level.results_fnames);

contrast=3;
cluster_forming_threshold=2.8;
num_perms=100;
corrp = osl_clustertf(permute(stats.lower_level_copes{contrast},[2 4 3 1]),cluster_forming_threshold,num_perms,26);

figure;
con=3; imagesc(stats.times, stats.frequencies, squeeze(corrp));axis xy;
ylabel('frequency (Hz)'); xlabel('time (s)'); colorbar; title(['cope' num2str(con)]);

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

disp('*************************************************');
disp('Finished osl_test_scrip_group test');
disp('*************************************************');
